home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / FORTRAN1.LZH / GAUSS.FOR < prev    next >
Text File  |  1988-02-08  |  7KB  |  265 lines
  1.       SUBROUTINE GAUSS ( A, Y, COEF, N, ERROR )
  2. C*
  3. C*                  *******************************
  4. C*                  *******************************
  5. C*                  **                           **
  6. C*                  **          GAUSS            **
  7. C*                  **                           **
  8. C*                  *******************************
  9. C*                  *******************************
  10. C*
  11. C*     SUBPROGRAM :
  12. C*          GAUSSIAN ELIMINATION
  13. C*
  14. C*     AUTHOR :
  15. C*          ART RAGOSTA
  16. C*          MS 207-5
  17. C*          AMES RESEARCH CENTER
  18. C*          MOFFETT FIELD,  CALIF  94035
  19. C*          (415)694-5578
  20. C*
  21. C*     PURPOSE :
  22. C*          TO SOLVE A SET OF SIMULTANEOUS, LINEAR EQUATIONS.
  23. C*          THE FORM OF THE EQUATIONS IS :
  24. C*                 Y1  =  A1,1*X1 +  A1,2*X2  + ...   A1,N*XN
  25. C*                 Y2  =  A2,1*X1...
  26. C*                 .
  27. C*                 .
  28. C*                 YN  =  AN,1*X1...                  AN,N*XN
  29. C*
  30. C*          THE SOLUTION IS OF THE FORM :
  31. C*                 X1 = COEF(1)
  32. C*                 X2 = COEF(2)
  33. C*                 .
  34. C*                 .
  35. C*                 XN = COEF(N)
  36. C*
  37. C*          REFERENCE : "PASCAL PROGRAMS FOR ENGINEERS AND SCIENTISTS",
  38. C*                      BY ALAN R. MILLER, 1981, SYBEX INC.
  39. C*
  40. C*     INPUT ARGUMENTS :
  41. C*          A     = N*N INPUT MATRIX(DESTROYED)
  42. C*          Y     = INPUT VECTOR OF LENGTH, N
  43. C*          N     = NUMBER OF EQUATIONS
  44. C*
  45. C*     OUTPUT ARGUMENTS :
  46. C*          COEF  = SOLUTION VECTOR OF LENGTH, N
  47. C*          ERROR = BOOLEAN ERROR FLAG (MATRIX SINGULAR IF TRUE)
  48. C*
  49. C*     INTERNAL WORK AREAS :
  50. C*          NONE
  51. C*
  52. C*     COMMON BLOCKS :
  53. C*          NONE
  54. C*
  55. C*     FILE REFERENCES :
  56. C*          NONE
  57. C*
  58. C*     SUBPROGRAM REFERENCES :
  59. C*          NONE
  60. C*
  61. C*     ERROR PROCESSING :
  62. C*          IF A ZERO APPEARS ON THE DIAGONAL AND CAN'T BE REMOVED, THE
  63. C*          MATRIX IS SINGULAR.
  64. C*
  65. C*     TRANSPORTABILITY LIMITATIONS :
  66. C*          NONE
  67. C*
  68. C*     ASSUMPTIONS AND RESTRICTIONS :
  69. C*          NONE
  70. C*
  71. C*     LANGUAGE AND COMPILER :
  72. C*          ANSI FORTRAN 77
  73. C*
  74. C*     VERSION AND DATE :
  75. C*          VERSION I.0      5-MAR-85
  76. C*
  77. C*     CHANGE HISTORY :
  78. C*           5-MAR-85    INITIAL VERSION
  79. C*
  80. C***********************************************************************
  81. C*
  82.       DIMENSION A(N,N), Y(N), COEF(N)
  83.       LOGICAL ERROR
  84. C
  85.       ERROR = .FALSE.
  86.       N1    = N - 1
  87.       DO 50 I = 1, N1
  88.          AMAX=ABS(A(I,I))
  89.          L = I
  90.          I1 = I+1
  91. C
  92. C ----- FIND LARGEST ELEMENT IN THIS COLUMN
  93. C
  94.          DO 10 J = I1, N
  95.             IF (ABS(A(J,I)) .GT. AMAX)THEN
  96.                AMAX = ABS(A(J,I))
  97.                L = J
  98.             ENDIF
  99. 10          CONTINUE
  100. C
  101. C ----- IF THE LARGEST ELEMENT IS ZERO, THE MATRIX IS SINGULAR
  102. C
  103.          IF (AMAX .EQ. 0.0)THEN
  104.             ERROR = .TRUE.
  105.             RETURN
  106.          ELSE
  107. C
  108. C -------- IF THE LARGEST ELEMENT IS NOT ALREADY ON THE DIAGONAL, PUT IT
  109. C --------  THERE BY SWAPPING ROWS
  110. C
  111.             IF (L .NE. I) THEN
  112.                DO 20 J = 1, N
  113.                   TEMP = A(L,J)
  114.                   A(L,J) = A(I,J)
  115.                   A(I,J) = TEMP
  116. 20                CONTINUE
  117.                TEMP = Y(L)
  118.                Y(L) = Y(I)
  119.                Y(I) = TEMP
  120.             ENDIF
  121. C
  122. C -------- DIVIDE EACH ELEMENT IN ROW BY THE LARGEST
  123. C
  124.             DO 40 J = I1, N
  125.                TEMP = A(J,I)/A(I,I)
  126. C
  127. C -------- NOW SUBTRACT THIS ROW FROM EACH SUBSEQUENT ROW
  128. C
  129.                DO 30 K = I1, N
  130.                   A(J,K) = A(J,K) - TEMP*A(I,K)
  131. 30                CONTINUE
  132.                Y(J) = Y(J) - TEMP*Y(I)
  133. 40             CONTINUE
  134.          ENDIF
  135. 50       CONTINUE
  136. C
  137. C --- IF A ZERO IS LEFT ON THE DIAGONAL, IT'S SINGULAR
  138. C
  139.       IF (A(N,N) .EQ. 0.0)THEN
  140.          ERROR = .TRUE.
  141. C
  142. C --- SUBSTITUTE FOR SOLUTION
  143. C
  144.       ELSE
  145.          COEF(N) = Y(N)/A(N,N)
  146.          I = N1
  147. 60       SUM = 0.0
  148.          I1 = I + 1
  149.          DO 70 J = I1, N
  150.             SUM = SUM + A(I,J)*COEF(J)
  151. 70          CONTINUE
  152.          COEF(I) = (Y(I)-SUM)/A(I,I)
  153.          I = I - 1
  154.          IF (I .GT. 0)GO TO 60
  155.       ENDIF
  156.       RETURN
  157.       END
  158. C
  159. C---END GAUSS
  160. C
  161.       CHARACTER FUNCTION GETC ( NIN )
  162. C*
  163. C*                  *******************************
  164. C*                  *******************************
  165. C*                  **                           **
  166. C*                  **          GETC             **
  167. C*                  **                           **
  168. C*                  *******************************
  169. C*                  *******************************
  170. C*
  171. C*     SUBPROGRAM :
  172. C*          GET CHARACTER
  173. C*
  174. C*     AUTHOR :
  175. C*          ART RAGOSTA
  176. C*          MS 207-5
  177. C*          AMES RESEARCH CENTER
  178. C*          MOFFETT FIELD, CA  94035
  179. C*          (415) 694-5578
  180. C*
  181. C*     PURPOSE :
  182. C*          GET A SINGLE CHARACTER FROM THE INPUT UNIT...
  183. C*            TAKE CARE OF NEW LINES AND END-OF-FILE
  184. C*
  185. C*     INPUT ARGUMENTS :
  186. C*          NIN - INPUT UNIT NUMBER
  187. C*
  188. C*     OUTPUT ARGUMENTS :
  189. C*          GETC - THE NEXT CHARACTER (FUNCTION VALUE)
  190. C*                 NOTE: CHAR(26) IS RETURNED FOR ENDFILE, CHAR(13) FOR ENDLINE
  191. C*
  192. C*     INTERNAL WORK AREAS :
  193. C*          NONE
  194. C*
  195. C*     COMMON BLOCKS :
  196. C*          NONE
  197. C*
  198. C*     FILE REFERENCES :
  199. C*          NIN
  200. C*
  201. C*     SUBPROGRAM REFERENCES :
  202. C*          MERROR, LENGTH
  203. C*
  204. C*     ERROR PROCESSING :
  205. C*          NONE
  206. C*
  207. C*     TRANSPORTABILITY LIMITATIONS :
  208. C*          THE SAVE STATEMENT (WHICH IS COMMENTED) MAY HAVE TO BE INCLUDED
  209. C*           ON SOME SYSTEMS.
  210. C*
  211. C*     ASSUMPTIONS AND RESTRICTIONS :
  212. C*          ALL TRAILING BLANKS ARE STRIPPED OFF
  213. C*          INPUT LINE LENGTH IS RESTRICTED TO 133 CHARACTERS
  214. C*
  215. C*     LANGUAGE AND COMPILER :
  216. C*          ANSI FORTRAN 77
  217. C*
  218. C*     VERSION AND DATE :
  219. C*          VERSION I.0     10-SEP-85
  220. C*
  221. C*     CHANGE HISTORY :
  222. C*          10-SEP-85    INITIAL VERSION
  223. C*
  224. C***********************************************************************
  225. C*
  226.       PARAMETER (LMAX=133,NSTART=LMAX+1)
  227.       CHARACTER *(LMAX) LINE
  228.       LOGICAL EOF
  229. C     SAVE NPTR, LPTR, EOF, LINE
  230.       DATA NPTR/NSTART/, LPTR/LMAX/
  231.       DATA EOF/.FALSE./
  232. C
  233. C --- END OF LINE WAS REACHED ON LAST ENTRY... GET NEW LINE
  234. C
  235.       IF (NPTR .GT. LPTR) THEN
  236.          IF ( EOF ) THEN
  237.             CALL MERROR('Attempted get after end of file.')
  238.             STOP
  239.          ENDIF
  240.          READ(NIN,900,END=1000) LINE
  241.          LPTR = LENGTH(LINE)+1
  242.          NPTR = 1
  243.       ENDIF
  244. C
  245. C --- RETURN NEXT CHARACTER UNLESS THERE ARE NO MORE... THEN RETURN <CR>
  246. C
  247.       IF (NPTR .EQ. LPTR) THEN
  248.          GETC = CHAR(13)
  249.       ELSE
  250.          GETC = LINE(NPTR:NPTR)
  251.       ENDIF
  252.       NPTR = NPTR + 1
  253.       RETURN
  254. C
  255. C --- END OF FILE
  256. C
  257. 1000  GETC = CHAR(26)
  258.       EOF = .TRUE.
  259.       RETURN
  260. 900   FORMAT(A)
  261.       END
  262. C
  263. C---END GETC
  264. C
  265.